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Abstract 

We discuss the crystallization process from the supersaturated melt in terms of its non¬ 
equilibrium properties. In particular, we quantify the amount of heat that is produced irreversibly 
when a suspension of hard spheres crystallizes. This amount of heat can be interpreted as arising 
from the resistance of the system against undergoing phase transition. We identify an intrinsic 
compression rate that separates a quasi-static regime from a regime of rapid crystallization. In the 
former the disspated heat grows linearly in the compression rate. In the latter the system crys¬ 
tallizes more easily, because new relaxation channels are opened, at the cost of forming a higher 
fraction of non-equilibrium crystal structures. In analogy to a shear-thinning fluid, the system 
shows a decreased resistance when it is driven rapidly. 
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Crystallization from the metastable melt is a non-equilibrium process. Yet, it is usually 
discussed in terms of quasi-equilibrium concepts such as transition state theorjff^^, which do 
not account for the fact that any irreversible process of hnite duration is inevitably subject 
to dissipation. Here, we present a numerical approach to assess the amount of energy that is 
dissipated in a crystallization process, and we discuss the relation between dissipation and 
external driving. 

The most obvious way to characterize the irreversibility of a process is to quantify entropy 
production. However, to do so directly is unpractical even for very simple model systems. 
To bypass this problem, we instead evaluate the mechanical work performed on the system 
by compressing it at a constant rate, and subtract the equilibrium work, which we obtain 
independently via the equation of state. 

As a model system we use hard spheres, the most simple system that shows a liquid- 
to-crystal transitioiP. Despite their simplicity hard spheres capture the essential physics 
of the crystallization process in many atomic systems. In general, at high densities the 
excluded volume between atoms dominates their dynamical behaviour, because the typical 
interparticle distances are smaller than the attraction ranges. The attractive forces thus 
effectively only add up to a flat background that does not influence the dynamics, but 
merely changes the equilibrium equation of statd^. Thus the results from our study should 
be applicable to a large class of crystallization processes in metallic systems as well as 
colloids. 

To model the crystallization process, we perform computer simulations of hard spheres in 
the NPT ensemble. The system is prepared in a fluid equilibrium state at constant pressure 
and then subjected to an increase in pressure with a constant rate P for a duration r. Under 
these conditions the work W performed on the system is 

W= r dtPv(t), 

Jo 

where V{t) is the volume response of the system to the external driving P. Assuming the 
difference in Gibbs free energy AG between the initial and the hnal state to be known, the 
dissipated energy of each simulation trajectory is 

IUdiss = fU-AG. 

The volume evolution of a typical trajectory is shown as a solid black line in Fig. (the 
volume is given in units of the sphere diameter cubed, a^, time in units of the free diffusion 
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FIG. 1. Evolution of the specific volume along a typical simulation trajectory of a hard-sphere 
system that is compressed with P = 0.01065 ksT/a^tQ (solid lines). Induction times are labeled 
t]\f and t'pf. Dashed lines indicate the equations of state of the fluid (upper line) and the ideal 
equilibrium crystal (lower line), respectively. A vertical line marks the time when the coexistence 
pressure Pc is crossed. Shaded areas indicate the different contributions to the work performed on 
the system, as discussed in the main text. Area (II) is the dissipated heat q^. 


time to which we specify in the next paragraph). Dashed lines indicate the equations of 
state of the liquicP and the crystaP. The dissipated energy ITdiss is indicated by the shaded 
areas. It consists of three contributions: (i) is the work associated with compression of the 
metastable liquid phase until the nucleation event occurs; (hi) is a contribution that arises 
because the system is not completely transformed into the equilibrium crystal during the 
simulation time, but contains defects. This contribution results in an almost constant offset 
in volume with respect to the equation of state of the equilibrium crystal, which does not 
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vary much between trajectories for any given P. Contribution (ii) yields the irreversible 
heat per particle, q^, associated with the crystallization process, 

rtN+At 

We model the dynamics of the system as stochastic, thus q'^ is a fluctuating quantity on the 
ensemble of trajectories. 

Before we discuss the distribution of we summarize the technical details of the 
simulation. We perform standard NPT Monte Carlo (MC) simulations with small par¬ 
ticle displacements drawn from a flat distribution from the interval [—A, A] with A = 
0.065(7. As unit of time we use to = cr^Z-Do, where the free-particle diffusion coefficient 
is Dq = A^/6/MC sweep ps 7 x 10“^(7^/MC sweep. To control the pressure, a volume 
change is attempted once per MC sweep by rescaling the box lengths according to Lj i—>■ 
Li exp[0.0012(r — 1/2)] where r is a uniform random variable in ]0,1] and i labels the Carte¬ 
sian directions. We allow changes of L* independently in each direction to accomodate 
crystals with unit cells of different aspect ratios. Simulations start from an equilibrium fluid 
state at pressure Pq = SksT/a^ and end at Pr = 23/cbT/( 7^, where ksT is the thermal 
energy and r = (P^ — Po)/|P| is the duration of the trajectory. The crystal-liquid coexis¬ 
tence pressure is Pc = 11.54 /cbT/( 7®. We monitor the degree of crystallinity by means of 
the local q^q^ bond order parametei®^. In order to distinguish different crystal structures 
we analyze the averaged bond order parameters \qi\ and Igep*- 

To compute q'^ we need to dehne the time window of contribution (ii) (see Fig. [^. We set 
the induction time tjsi to the time after which the largest crystalline cluster maintains a size 
of ten or more particles. The end of the process, tvf + At, is set to the time when the overall 
crystallinity reaches 60%. This value is large enough to capture the main contributions 
of dissipated heat, but still small enough to minimize the influence of periodic boundary 
conditions. 

Since rare trajectories can contribute considerably to the non-equilibrium work distribu¬ 
tion, we need to generate a very large number of independent trajectories. As the com¬ 
putational effort is large (0(10®) trajectories per value of compression rate), we simulate a 
relatively small system of A = 540 particles, a system size that, albeit small, still repro¬ 
duces the nucleation rates correctly. To verify this, we compared the nucleation rates to 
those obtained from simulations with N = 8000 and N = 216, 000 particle^. 
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We compressed the system for times r = 1 x 10^, 2 x 10^, 5 x 10®, 1 x 10®, 2 x 10®, 5 x 10®, 
and 1 X 10^ MC sweeps (corresponding to compression rates between P ^ 0.214 ksT/aHo 
and P = 0.00214: ksT/aHo.) The number of trajectories sampled varied between 70,000 
and 650, 000 depending on the accuracy needed for a given value of P. In total this required 
90 years of CPU time on 2.2GHz Xeond^. 

The left panel of Fig. shows the distribution of work performed on the ensemble of 
non-equilibrium trajectories. Solid lines mark p{W/N) for different values of |P| in the 
forward process {pf{W/N) on compression, P > 0) and the reverse process (pniW/N) on 
expansion shown as pui—W/N), P < 0,). The distributions for the expansion processes are 
centered around values |IU|/iV < |A/i|, where A/r = AG/N is the difference in chemical 
potential between the initial and the hnal state. For all values of P < 0, the curves are well 
described by Gaussian probability distributions down to the accuracy set by the number of 
trajectories that we simulated. The distributions pf{W/N) associated to the compression 
processes are centered around W/N > A/i. In particular at small P, they deviate from 
Gaussian behaviour and display a more subtle structure which we discuss later in terms of 
q^. 

To estimate whether our ensemble averaging is sufficient, we test whether these work 
distributions are compatible with the Jarzinsky relationl^l and the Grooks relation!!^. For 
an arbitrary non-equilibrium process, the Jarzinsky relation connects the work distribution 
p{W) to the equilibrium free energy difference AG, 

{exp{-/3W)) = jp{W) exp(-/3IU) dW = exp(-/3AG). (1) 

Since the average is over the exponential of IF, the Jarzynski relation provides a very 
sensitive test for the accuracy of sampling. 

For distributions p{W) that are superpositions of Gaussian distributions, eq. can be 
evaluated analytical!}^. Fitting the work distribution of the forward process with a super¬ 
position of two Gaussians, constrained such that the backward process is also well described 
by Grooks’ fluctuation theorem, we use the distributions p{W) from our simulations to es¬ 
timate A/r. The results are shown in the right panel of Fig. (circles with error bars). 
Agreement with the equilibrium chemical potential difference obtained from the equation of 
state is reasonable, thus we conclude that our sampling is sufficient. 

As discussed above, most of the work IF performed on the system during compression 
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FIG. 2. Left panel: Distribution of work per particle, p{W/N), performed upon compression 
(forward process: right set of curves) and expansion (reverse process: left set of curves, shown 
as p{—W/N)) across the phase transition with a constant rate |P|. Histograms are shown for 
|P| = 0.214,0.107,0.0428,0.0107,0.00428/cT/cj^to (right to left for forward process). The dashed 
vertical line indicates the equilibrium chemical potential difference A|U. Right panel: A/r estimated 
using the Jarzynski relation, eq. see text. The horizontal dashed line indicates the equilibrium 
value. 

consists of equilibrium or quasi-equilibrium contributions that are readily evaluated if one 
knows the equations of state of the initial and dual phase. The non-equilibrium nature of 
the process is characterized by the distribution of dissipated energy; a quantity that has not 
been discussed before in the literature on crystallization. 

The left panel of Fig. shows the distribution of dissipated energy per particle, for 
various values of the compression rate P. As the compression rate increases, the distribution 
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FIG. 3. Left panel: Probability distribution of dissipated energy per particle during crystallization, 
q'^, for different compression rates, P = 0.00214, 0.0107,0.0214, 0.107, 0.214/cT/cr^fO) as indicated. 
Right panel: Corresponding distributions of the crystallization loss, q^/P. 

shifts to higher average {q'^) and broadens. At the highest compression rate we simnlated, 
{qp is abont 0.2 which is of the same order of magnitnde as the average (macroscopic) 
interfacial energy over the area per particle, ~ 0.6 In the right panel of the hgure 

we show that the distribntions collapse for weak driving P, when plotted in terms of the 
reduced variable g'^/ P. This collapse defines the regime of qnasi-static behavionr, where the 
response q^/P of the system is independent of the driving force. 

In the context of eqnilibrinm thermodynamics the term “qnasi-static” is restricted to 
the case of inhnitely slow driving, P = 0. The existence of a regime of P-independent 
distribntions of g'^/P jnstifies the extension of this notion to finite (small) driving rates. The 
limiting value of the average response C, := {qp/P attained for P —)■ 0 can be interpreted 
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FIG. 4. Crystallization loss, i.e,, the average energy per particle dissipated during crystallization 
relative to the external driving rate,(g'^)/P as a function of compression rate P. 

as an immanent system property, the quasi-static crystallization loss. For the hard-sphere 
system, we obtain Cp^o ~ l-6cr^to- 

At driving forces above a threshold P* the crystallization loss C{P), as a fnnction of P, 
drops sharply (see Fig. |^.(If one took into acconnt the energy costs due to defects in the 
crystal (contribution (iii) in Fig. this effect would even be enhanced, since the excess 
volume over the equilibrium volume Kq (used to dehne q'^) increases monotonically with 
increasing P.) Intuitively, one would expect the relative dissipation to increase once the 
rate of driving exceeds typical microscopic relaxation times of the system, as additional 
work can be dissipated through the microscopic degrees of freedom. The counter-intuitive 
behaviour of the crystallization loss can be rationalized by analogy with mechanical friction 
in fluids. There, one typically observes friction to decrease strongly in the nonlinear-response 
regime of fast drivin^^^. This is particularly well known for the viscosity of non-Newtonian 
fluid^, where the effect is called shear thinning. It also holds for a driven tracer subject 
to an external force in a dense fluicP®. In these cases, the slow near-equilibrium relaxation 
processes are replaced by faster ones that occur on the time scale set by the external driving. 
In analogy, we interpret C{P) as a generalized friction coefficient that characterizes the melt’s 
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resistance to phase transformation. 

C{P) shows non-monotonic behavionr becanse both effects, i.e., increased friction throngh 
enhanced conpling to microscopic degrees of freedom as well as decreased friction throngh 
non-eqnilibrinm relaxation channels, contribnte to the crystallization loss. As indicated in 
Fig. the initial increase is associated with an increase in the large-g'^ tail. This is intnitively 
expected since an increased conpling to microscopic degrees of freedom increases the proba¬ 
bility for stronlgy dissipating trajectories. At P > P*, this large-g^ tail is cnt off. (We will 
show in the following that this effect is dne to the formation of non-eqnilibrinm crystal strnc- 
tnres.) The crossover between the two trends occnrs aronnd P* ^ 2 x 10“^ fcsT/cr^to-This 
cross-over valne is explained by the time scale needed for collective particle rearrange¬ 
ments involving the nearest and next-to-nearest neighbonr shells, ti is set by the long-time 
self-diffnsion coefficient Dl = jti^. For the typical densities reached when crystallization 
sets in, Dl/Dq = 0(10“^) (for the initial flnid state in onr work, Dl/Dq ^ 0 04)1211. Hence, 
P*tL = OlksT/cr^)] i.e. the effects of the external driving start to dominate the crystalliza¬ 
tion process once the compression rate is faster than the typical thermal energy density can 
be redistribnted throngh collective particle rearrangements. 

To demonstrate that the melt indeed relaxes faster into the crystal phase at P > P*, 
we show in Fig. |^the distribntions of the crystallization time At (i.e. the distribntions of 
the length of time between the indnction time and the time when 60% of the system are 
crystallized). For small P, the distribntions again collapse to a P-independent cnrve. This 
curve displays a pronounced tail at large At, and the crystallization process is slow on the 
time scale to of free particle diffusion. The average (At) is approximately 20to, i.e. on the 
order of the long-time self-diffusion time t^. This fact conhrms that long-time diffusion sets 
the relevant time scale for the crystallization process. At large P, the distributions shift to 
smaller average At, and the large-At tail is cut off. Moreover, p(At) narrows proportionally 
to P, which suggests that the inverse external driving rate 1/P sets the relevant time scale 
for the dynamics. The change in shape of the distribution p{At) with P is qualitatively 
similar to the one observed for the distribution p{q'^/P) shown in Fig. This emphasizes 
that the change in dissipation mechanism is of kinetic rather than thermodynamic origin. 

Next we show that the accelerated crystallization mechanism proceeds through non- 
equilibrium relaxation channels, in particular the formation of non-equilibrium crystal struc¬ 
tures (bcc instead of fee). Fig. [^shows the probability distribution of the fraction of particles 
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FIG. 5. Distribution of the time interval over which the system crystallizes, for different P 
corresponding to the data shown in Fig. Vertical lines indicate the maximum At possible in the 
simulation for the earliest induction time observed in our simulations, for the highest two P 
shown. 

with a bcc-like environment in the crystal at the end of the simulation run. Again, at slow 
driving rates, P < P*, the distributions are independent of P. For P > P*, more bcc 
structures are formed. There even is a signihant number of runs that crystallize completely 
into bcc. Our data indicate that the formation of bcc-like structures is facilitated at large 
compression rates. Interestingly, the question why fee is the stable equilibrium structure 
while bcc should form more easily is known from Landau theory®. There, the effect arises 
because a larger set of reciprocal lattice vectors is needed to form fee. This implies that a 
larger set of local density fluctuations needs to be sampled. It is conceivable that this takes 
more time, and hence, bcc is favored kinetically. 

The tendency to form metastable crystal structures in rapid solidiheation is well known 
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FIG. 6. Distribution of the fraction of bcc crystal structures in the crystalline part of the system 
at the end of the compression run, for different P corresponding to the data shown in Fig. 

from metallic melt^. It is often attributed to Oswald’s step rule, which invokes interfaction 
tensions between the crystal nucleus and the surrounding fluid. We offer an alternative 
explanation that is founded on microscopic kinetic arguments, rather than macroscopic 
thermodynamic quantities that might not be well defined on the scale of a few particle 
diameters. 

CONCLUSION 

We have discussed crystallization in terms of non-equilibrium notions and calculated the 
distribution of heat dissipated during a crystallization process. Compressing the system at 
different rates P, we measure the volume response and hnd two regimes: Below a charac¬ 
teristic compression rate, P*, set by the single-particle diffusion time and the coexistence 
pressure, the resistance of the system against the phase transition, {qp/P, is constant. The 
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system responds quasi-statically. Above P* the crystallization process evolves far from equi¬ 
librium. The system crystallizes more easily than expected, because new relaxation channels 
are opened via the formation of bcc structures instead of the thermodynamically favored 
fee ones. In this regime the evolution of the system is determined by kinetics rather than 
thermodynamics. 
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